'''
Created on Sep 8, 2010

@author: oabalbin
'''

import os
import sys
import glob
import subprocess
import pysam


import numpy as np
#import matplotlib.pyplot as plt
'''
import matplotlib.mlab as mlab
import matplotlib.pylab as pla
import matplotlib.cm as cm 
'''

#import rpy2.robjects as robjects
#import rpy2.robjects.numpy2ri
#from rpy2.robjects.packages import importr

from Bio.Seq import Seq
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import generic_dna


def extract_sequences_from_samfile(fastq_file_name):
    
    mysam_file = pysam.Samfile(fastq_file_name,'r')

    reads_quality_scores, reads_trimming_length=[],[]
    reads_quality_scores_str, reads_trimming_length_str=[],[]
    mybins=np.arange(0,41,1)
    for alignedread in mysam_file.fetch():         

        myletter_annotations={}
        # Check the header for TopHat @PG    ID:TopHat    VN:1.0.13
        if mysam_file.header['PG']['ID']=='TopHat' and mysam_file.header['PG']['VN']=='1.0.13':
            # Note that the sam file says to store the base quality score as ASCII-33.
            myletter_annotations['phred_quality'] = [ord(letter)-33 for letter in list(alignedread.qual)]        
        else:
            # Need to see what to do with other versions.
            print fastq_file_name, mysam_file.header['PG']
            continue
        
        rec = SeqRecord(Seq(",".join(list(alignedread.seq)).replace(',',''), generic_dna), id=alignedread.qname, name=alignedread.qname, description="sequence from a sam file",\
                        letter_annotations=myletter_annotations)
        rec.format="fastq"
    
        print rec.letter_annotations["phred_quality"]
        print alignedread.cigar
        print alignedread
    #    if c < 100000:
    #    print rec.id
        trim=len(rec.letter_annotations["phred_quality"])
        # Check if there is something similar to B (fastq-illumina) in the Phred Sanger format
        '''
        try:
            trim = rec.letter_annotations["phred_quality"].index(2)
        except ValueError:
            pass
        '''
        reads_quality_scores.append(np.mean(rec.letter_annotations["phred_quality"][:trim]))
        reads_trimming_length.append(len(rec.letter_annotations["phred_quality"][:trim]))
        ##
        reads_quality_scores_str.append(str(np.mean(rec.letter_annotations["phred_quality"][:trim])))
        reads_trimming_length_str.append(str(len(rec.letter_annotations["phred_quality"][:trim])))

     #       c+=1
    
    mysam_file.close()
    
    return reads_quality_scores, reads_trimming_length, reads_quality_scores_str,reads_trimming_length_str



fastq_file_name = '/data/projects/mirnas/mirnas_discovery/sam_files/test/mctp_30Y5NAAXX_6_test_tu_complilation.sam'#RB1.sam'
print fastq_file_name

extract_sequences_from_samfile(fastq_file_name)